%This code produces three plots presenting the comparative statics wrt phi

%Last update 03/01/2023
clear all
close all
clc


myopt = optimset('Display','Off','MaxFunEvals',1e1000000,'TolFun',1e-20,'TolX',1e-30);

datadir = '..\..\Data\Final\Figures11_12\';
figdir = '..\..\Plots\';

warning off
load([datadir 'phi_data_S100J50_alpha_0001_ADJalpha']);

%%

vec_q=[grid_q(1) 1 grid_q(end)];
q3=vec_q(3);
q2=vec_q(2);
q1=vec_q(1);

for qq=1:length(vec_q)
    q=vec_q(qq);
    for m=1:M
        for s=1:S
            n_s = pchip(grid_q,n_mat_save(:,s,m));
            sales_s = pchip(grid_q,sales_mat_save(:,s,m));
            plants(qq,s,m) = ppval(n_s,q);
            sales(qq,s,m) = ppval(sales_s,q);
        end
    end
end

fig1 = figure(1);
qh=plot(log(grid_income),plants(3,:,1),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),plants(2,:,1),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),plants(1,:,1),'linewidth',2.5,'Color',[91, 207, 244] / 255);
hold on
qh=plot(log(grid_income),plants(3,:,M),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
hold on
qh=plot(log(grid_income),plants(2,:,M),'linewidth',2.5,'Color',[0 0 1],'linestyle',':');
hold on
qh=plot(log(grid_income),plants(1,:,M),'linewidth',2.5,'Color',[91, 207, 244] / 255,'linestyle',':');
th=legend("$q=$"+q3+", $\phi=$"+run_vec(1)+"","$q=$"+q2+", $\phi=$"+run_vec(1)+"","$q=$"+q1+", $\phi=$"+run_vec(1)+"","$q=$"+q3+", $\phi=$"+run_vec(M)+"","$q=$"+q2+", $\phi=$"+run_vec(M)+"","$q=$"+q1+", $\phi=$"+run_vec(M)+"",'Location','best');
set(th,'Fontsize',9,'Interpreter','Latex');
th=ylabel("Measure of Plants of Firm $j$ in $s$, $n_{js}$");
set(th,'Fontsize',12,'Interpreter','Latex');
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex');
axis tight
saveas(fig1, [figdir 'Figure11_left'], 'png');
%%

for m=1:M
    for j=1:length(vec_q)
        for s=1:S
            if sales(j,s,m)<=0
                sales_plot(j,s,m)=0.00000001;
            else
                sales_plot(j,s,m)=sales(j,s,m);
            end
        end
    end
end

fig2 = figure(2);
qh=plot(log(grid_income),log(sales_plot(3,:,1)),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_income),log(sales_plot(2,:,1)),'linewidth',2.5,'Color',[0 0 1]);
hold on
qh=plot(log(grid_income),log(sales_plot(1,:,1)),'linewidth',2.5,'Color',[91, 207, 244] / 255);
hold on
qh=plot(log(grid_income),log(sales_plot(3,:,end)),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
hold on
qh=plot(log(grid_income),log(sales_plot(2,:,end)),'linewidth',2.5,'Color',[0 0 1],'linestyle',':');
hold on
qh=plot(log(grid_income),log(sales_plot(1,:,end)),'linewidth',2.5,'Color',[91, 207, 244] / 255,'linestyle',':');
th=legend("$q=$"+q3+", $\phi=$"+run_vec(1)+"","$q=$"+q2+", $\phi=$"+run_vec(1)+"","$q=$"+q1+", $\phi=$"+run_vec(1)+"","$q=$"+q3+", $\phi=$"+run_vec(M)+"","$q=$"+q2+", $\phi=$"+run_vec(M)+"","$q=$"+q1+", $\phi=$"+run_vec(M)+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex');
th=ylabel("Log of Sales in $s$");
set(th,'Fontsize',12,'Interpreter','Latex');
th=xlabel("Log of Income in $s$, $\ln I_s$");
set(th,'Fontsize',12,'Interpreter','Latex')
axis([log(grid_income(1)) log(grid_income(end)) -2 max(log(sales_plot(end,:,end)))]);
saveas(fig2,[figdir 'Figure11_right'],'png');


%% Markets

for j=1:J
    if Total_sales_j(j,end)==0
        tot_sales_plot_end(j,1)=0.000001;
    else
        tot_sales_plot_end(j,1)=Total_sales_j(j,end);
    end
end

fig3 = figure(3);
qh=plot(log(grid_q),log(Total_sales_j(:,1)),'linewidth',2.5,'Color',[0 0 0]);
hold on
qh=plot(log(grid_q),log(tot_sales_plot_end(:,end)),'linewidth',2.5,'Color',[0 0 0],'linestyle',':');
th=legend("$\phi=$"+run_vec(1)+"","$\phi=$"+run_vec(M)+"",'Location','northwest');
set(th,'Fontsize',9,'Interpreter','Latex');
th=ylabel("Log of Total Sales of Firm $j$");
set(th,'Fontsize',12,'Interpreter','Latex');
th=xlabel("Log of Firm $j$ Productivity, $\ln q_j$");
set(th,'Fontsize',12,'Interpreter','Latex');
axis([log(grid_q(1)) log(grid_q(end)) -2 log(max(Total_sales_j(end,1)',tot_sales_plot_end(end,end)'))]);
saveas(fig3, [figdir 'Figure12'],'png');


